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Abstract: In recent years, various applications regarding sparse continuous signal recovery 
such as source localization, radar imaging, communication channel estimation, etc., have 
been addressed from the perspective of compressive sensing (CS) theory. However, there 
are two major defects that need to be tackled when considering any practical utilization. 
The first issue is off-grid problem caused by the basis mismatch between arbitrary located 
unknowns and the pre-specified dictionary, which would make conventional CS reconstruction 
methods degrade considerably. The second important issue is the urgent demand for 
low-complexity algorithms, especially when faced with the requirement of real-time 
implementation. In this paper, to deal with these two problems, we have presented 
three fast and accurate sparse reconstruction algorithms, termed as HR-DCD, Hlog-DCD 
and H/p-DCD, which are based on homotopy, dichotomous coordinate descent (DCD) 
iterations and non-convex regularizations, by combining with the grid refinement technique. 
Experimental results are provided to demonstrate the effectiveness of the proposed algorithms 
and related analysis. 

Keywords: compressed sensing; sparse continuous signal recovery; off-grid problem; 
low-complexity reconstruction; non-convex regularization 
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1. Introduction 

Compressive sensing (CS) and sparse signal representation have received widespread attention and 
increasing interest over the past few years [1,2], which are motivated by the sparse nature of the real 
world data and the advantages of the CS theory. The applications of CS in numerous areas have been 
widely investigated in the literature, such as magnetic resonance imaging (MRI) [3], synthetic aperture 
radar (SAR) imaging [4], inverse synthetic aperture radar (ISAR) imaging [5], passive radar imaging [6], 
direction-of-arrival (DO A) estimation [7], conmiunication channel estimation [8], seismic signal 
processing [9], spectral estimation [10], image processing [11], speech enhancement [12], etc. Generally 
speaking, those disciplines explore the following linear signal model [13]: 



where y G c^^^ is a M X 1 complex vector, H G c^^^ is a M X complex measurement matrix, 
w G c^^^ represents the X 1 unknown complex signals of interest, and n G d^^^ denotes a M X 1 
complex noise vector. Suppose w is sparse or compressible in a known dictionary V G c^^^, i.e., 
w = Vx, where x G c^^^ is a ^-sparse vector, namely it can be approximated by its K largest 
coefficients or its coefficients satisfy a power decay law with K strongest coefficients. Therefore, 
linear measurements are obtained in CS as: 



where A = WV G c^^^ represents the sensing matrix. Although the above equation is usually ill-posed, 
the CS theory has shown that if A satisfies some certain conditions, we can construct x and w stably 
from highly undersampled measurements y [14]. 

1.1. Off-Grid Problem in CS-Based Methods 

In the CS processing procedure, the first necessary step is to design a dictionary through a discretization 
operation with the assumption that the elements of unknown x lie exactly on those pre-defined grids 
corresponding to A. Obviously, this is practically impossible since the continuous nature of the 
unknowns of x, such as the unknown directions, may not fall into the predefined angular grids in DOA, 
and the scatterers of the target may not locate exactly on the pre-discretized scene grids in radar 
imaging. Hence, once faced with a continuous signal, the off-grid problem in using the conventional 
sparse recovery techniques is inevitable, no matter how densely we grid x. Previous researches [15-17] 
have demonstrated that the traditional CS-based methods would be severely affected when the off-grid 
problem emerges, and [18] also claimed that the off-grid problem is one of the major constraints in 
popularizing CS-based methods in practice. It is worth noting that there may be other factors leading to 
basis mismatch, for example, in the radar imaging field, unsatisfactory system parameters (e.g., 
position errors of transceivers [19] and phase synchronization mismatch [20]) are also likely to 
degrade the performance of conventional CS-based methods from our previous researches, however 
this paper only focuses on the off-grid problem, and assumes the system errors are small enough. 

The solutions for off-grid problem have been broadly studied in previous literatures [21-30]. So far, 
main research topics include three types as follows: 



y = Hw + n 



(1) 



y = Ax + n 



(2) 
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(a) Direct method according to theories of Xampling and the finite rate of innovation (FRI). This 
scheme is first introduced by Michaeh et al. as a general framework within various solutions for 
analog signals [21]. The significant advantage of this method is its direct solution, without any 
pre-discretization of the continuous signal at first. Thus it is less sensitive to off-grid problem 
compared with traditional CS-based methods. However, this method is mainly based on spectral 
estimation algorithms (e.g., ESPIRT [22], MUSIC [6], Matrix Pencil [23]), which needs special 
signal expression in linear measurement model and its performance may suffer from low SNR 
and few snapshots [24]. 

(b) Off-grid CS method. This kind of method is under a first order Taylor expansion model, and 
thus is quite sensitively depended on the model's accuracy, and cannot give a thorough solution 
when higher order approximations are significant [23]. Furthermore, it is likely to involve 
highly-computational burden [25,26] for alternatively finding the sparse solution and estimating 
the off-grid error, especially in large scale problem. 

(c) Grid refinement approach. The idea of the grid refinement technique was firstly introduced by 
Malioutov et al. [27] to mitigate the effect of limiting estimates to a grid of spatial locations in 
source localization problem. Then it is generalized by Liu et al. [28] for locating two-dimensional 
multiple underwater acoustic sources. Recently, Guldogan et al. [29] have proposed a novel grid 
refinement algorithm to alleviate off-grid problem by using a particle swarm optimization (PSO) 
and orthogonal matching pursuit (OMP), which makes use of the PSO to perturb the location of 
each grid point. 

Since the "grid refinement approach" is iteratively refined to match with the desired resolution of 
the off-grid components by "coarse and fine grid partition" operations, and does not have the 
drawbacks of (a) and (b), this paper follows the idea of "grid refinement approach" which makes a 
coarse grid first instead of having an universally fine grid to reduce the complexity, then achieves fine 
grids around the peaks using more refinement levels. In this way, after a few iterations of refining 
process, it becomes fine enough that off-grid problem effect is negligible. 

1.2. Fast and Efficient Algorithms for Real-Time Implementation 

Furthermore, there exists another common challenge by utilizing the CS-based methods, i.e., we 
need fast and efficient algorithms for real-time system implementation, particularly for digital electronic 
circuits (e.g., ARM, FPGA, DSP) [30]. From previous researches, sparse recovery techniques can be 
roughly divided into two families, that are greedy methods (e.g., MP [31], OMP [32], GP [33]) and 
optimization based methods (e.g., l\ optimization [34], smoothed /o optimization [35], non-convex 
optimization [36]). Generally speaking, greedy methods have the advantages of lower complexity, 
faster speed, less storage requirement, and flexible implementation compared with optimization 
based methods, and are considered the most suitable candidates for hardware implementation. 
However, their performances are inferior to those of the optimization based methods, such as h 
norm minimization basis pursuit de-noising (BPDN) [34]. 
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As the coordinate descent (CD) search has an inherent property of being low complexity when 
signals are sparse, Zakharov et al. have successfully exploited dichotomous CD (DCD) iterations for 
solving LS [37], RLS [38] and MVDR beamforming problems [39] using FPGAs for real-time 
implementation. In addition, to solve the reweighted h minimization problem, they have developed a 
greedy algorithm called H/i-DCD [40], i.e., homotopy DCD method with reweighted h penalty, which 
has shown a high recovery performance with relatively low complexity. It turns out that the overall 
complexity of the algorithm is comparable to that of MP. Moreover, homotopy iterations result 
in high accuracy of the solution, even higher than that of the YALLl algorithm. Meanwhile, 
Liu and Zakharov et al. [28] have utilized the low complexity homotopy approach combined with CD 
search for locating underwater acoustic sources by solving the multi-frequency BPDN problem. The 
proposed method is evaluated by applying to simulated and real experimental data. As far as we know, 
previous studies [28,37^0] are mainly focusing on homotopy DCD with convex regularizations, while 
there are seldom researches linking to homotopy DCD with non-convex regularizations. 

1.3. Our Contribution 

From the above discussion, we have addressed two major issues {i.e., off-grid problem and efficient 
algorithms for real-time implementation) of applying CS to sparse continuous signal reconstruction. 
In this paper, motivated by the idea of H/l-DCD [40] and grid refinement approach [27], we 
present three fast and accurate sparse reconstruction algorithms {i.e., HR-DCD, Hlog-DCD and 
H/p-DCD) which are based on homotopy, dichotomous coordinate descent (DCD) iterations and 
non-convex regularization, combining with the grid refinement technique to deal with the aforementioned 
issues. The main contributions of this paper are as follows: (1) we formulate the sparse recovery 
problem by homotopy DCD method with three typical classes of non-convex penalties, which are 
proved to recover sparsity in a more efficient way than homotopy DCD method with convex penalties 
as shown in Zakharov's previous researches [28,37^0]; (2) further, the grid refinement technique [27] 
is utilized to combine with our algorithms to alleviate the effect of off-grid problem and reduce the 
complexity simultaneously; (3) experimental results of three representative applications (DOA, passive 
radar imaging, ISAR imaging) are carried out to verify the effectiveness of the proposed methods. 

The outline of this paper is as follows: in Section 2, two sparse recovery problems are discussed. 
Section 3 describes the proposed methods and related analysis. In Section 4, extensive experimental 
results are presented to verify our methods. Finally, Section 5 draws the conclusions. 

Throughout this paper we shall make use of the following notation: bold-case letters are reserved 
for vectors and matrices, e.g., x is a vector, A is a matrix; Elements of A and x are represented as A„,p 
and x„, respectively. / and f denote the support and its complement, respectively. Further, we represent 
by R^^^ the q-ih column of R; A/ a matrix obtained from A keeping only columns corresponding to /; 
X/ the subset of x that contains entries from x corresponding to /; ||x||p denotes the Ip norm of a vector, 
II All 2 is the spectral norm of the matrix A; {, ) denotes the inner product; (O^, {'Y, Re{-} denote the 
conjugate transpose, conjugate and real part operations, respectively. 
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2. Problem Formulation and Motivations 

This section briefly introduces the general sparse recovery framework, upon which we develop 
our algorithms. As stated before, various applications can be represented by the linear model as 
Equation (2) shows. 

When considering the practical applications, sparse recovery techniques require low-complexity 
algorithms just as some kinds of greedy methods which are suitable for real-time implementation. 
However, the performances of greedy methods are compared unfavorably with optimization based 
methods (e.g., using l\ norm, reweighted l\ norm, Ip norm, smoothed /o regularization). Recently, 
Zakharov et al have proposed the homotopy DCD method with convex regularizations, including l\ 
and reweighted l\ norms, which have shown a high recovery performance and relatively low complexity. 
Moreover, most of the operations in the algorithms are additions, therefore, they are very suited for the 
hardware implementation of real-time operating systems. Motivated by the idea of homotopy DCD 
method and the fact that non-convex regularization usually yields a sparser solution than any convex 
penalty for a given residual energy (see Figure 1 for example), this paper proposes a fast and accurate 
sparse signal reconstruction by homotopy DCD technique with non-convex regularizations. Three typical 
non-convex penalty functions [12,41] are considered, Le,, the first order rational function penalty, the 
logarithmic penalty and the Ip penalty, respectively. With the corresponding penalties, we have 
achieved the derivation of the HR-DCD, Hlog-DCD and H/^-DCD algorithms in the next section. 



Figure 1. Examples of different penalty functions. 
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Since x is usually distributed continuously in the corresponding space in many applications, 
the off-grid problem is likely to exist. For example, in the radar imaging field, the reflectivity centers 
of the scatterers are generally not located at exact on-grid spatial positions illustrated in Figure 2, 
which means the measurement matrix should cover the basis vectors corresponding to off-grid 
scatterers. As conventional CS -based methods do not consider the off-grid effect, therefore the related 
sparse reconstruction techniques suffer unacceptable degradation in image quality. This paper has 
utilized the grid refinement approach [27] shown in Figure 3 combined with the HR-DCD, Hlog-DCD 
and H/p-DCD algorithms to alleviate the impact of off-grid problem (see details in Section 3.2). 
Experimental results in Section 4 are provided to demonstrate the performance improvement of the 
proposed algorithms. 
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Figure 2. Off grid and on grid scatterers. 

X X X X X 



X Grid 



• Unknowns 



Figure 3. Illustration of grid refinement approach. 
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3. Homotopy DCD with Non-Convex Regularization Algorithms 

We consider minimization of the following cost function to solve the problem Equation (2): 

1 2 



(1) 



where the sparsity measurement function p(x) is generally separable, or p(x) = Hi 0 (^i)' where 
denotes the parameterized form and satisfies such conditions that p(x) is sparsity-preserving, and 
T > 0 is a regularization parameter. As mentioned in previous researches, we can use DCD iterations 
for minimizing the cost function, which is represented in a fixed point. The DCD algorithm is 
appealing for practical designs as it operates at the bit level, resulting in stable hardware 
implementations, which is shown in Table 1. Different from Zakharov's studies [28,37^0], this paper 
utilizes the non-convex regularizations instead of convex regularizations embedding with the DCD 
algorithm, and chooses three typical penalty functions as follows: 
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(a) The first order rational function penalty: 



(2) 



(b) The logarithmic penalty: 



(c) The Ip penalty: 



(fi{x,) = -log{l^a\x,\) (3) 



x,) = \x,\' (4) 



In order to reduce the complexity of the DCD algorithm, we develop a greedy algorithm that is 
based on homotopy method with respect to a set of the parameter r. Besides, the updates of our DCD 
algorithm are only done within the support instead of all elements. As the support is usually much 
smaller than A^, therefore, the complexity is further reduced. We consider the first order rational 
function penalty at first, and develop HR-DCD algorithm. Then Hlog-DCD algorithm and H/^-DCD 
algorithm can be similarly obtained. 

Table 1. DCD Algorithm. 

Step Equation 

Initialization: x = 0, r = y, R = A^A, m = 0, S = H 

1 for m = 1 to Mb, repeat: 

2 S = 5/2, a = [5, -5,j5, -jS] 

3 for / = 1 to A^, repeat: 

4 for A: = 1 to 4, repeat: 

5 a; = S^Rt,i/2 - Re {<(A«)''r} + t(0(x, + aj,) - 0(x,)) 

6 ifA;<0,do: 

8 end if 

9 end for 

10 end for 

1 1 end for 



The following two propositions define the rules for adding/removing elements into/from the support 
of HR-DCD algorithm. 

Let / be the support at some homotopy iteration, and / be its complementary set. Denote r = y - Ax, 
R = A^A and b = A^r. 

Proposition 1: Add the ^-th (t E r) element into the support / according to the rule: 
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Proposition 2: Remove the ^th (t G /) element from the support / according to the rule 

? = argimn-|x f , |-r — ^— — s.t. -Ixf +9i(xZ?,|-r — ^— — <0 



1 + - X 



1 ^11 
1 + - X 

2' 



(6) 



We now prove Proposition 1 . 

Proof: Let the ^th element = 0 be activated as = a = and the updated solution 

vector is denoted as x. The update of the cost function Equation (3) is then given by: 



AJ = j[x]-j[x] = haf R,,-^{ab^} + T^ 



1 ^1 I 
1 + - d 

2' ' 



(7) 



The cost function is reduced if A/ < 0 . For a fixed \a\ , AJ achieves a minimum value if 
arg{a) = arg{/?t), and in this case: 



1 I |2 I II I \^ 

AJ = -\a\ K-\a\\bA-\-T — — 



1 + - a 
2' ' 



Let: 



AJ' = 4\ l + -\a\ 



Ay 



a\ aR^^ \a\ + 2[R^ ^ - a\b^\)\a\^ 4[t - \b^\) 



= aRAa\- 



r, ,R.,-a\b,\Y fR^^^-a\b,\]\4{T-\b,\) 



\a +- 



aR, 



aR, 



(8) 



(9) 



For \a\ > 0, if |a| = - we have: 



(10) 



2 

If AJ' < 0, then 4a7?t - (T?^ ^ + a|&t|) < 0. Thus, there exists a value of a that results in AJ < 0, 
in reducing the cost function. It is seen from the last expression that if we want to add a new 
element to the solution vector, the index t of the element should correspond to the maximum of 
~^~^^T~^((^t,t + ^l^tl)^ ~ 4a7?t over t G In this case, we will obtain the largest decrement of 

the cost function. 

Thus, the value of a that for a fixed t such that AaR^iT — {R^^ + a\bi\) < 0 results in the 
decrement of the cost function is given by: 



a\bA-R, 



aR, 



(11) 



According to the above statements, the (n + l)-th nonzero element to be activated in x should satisfy: 



t = ars 



For I a I > 0, we have: 
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1 I |2 I II I m 

A/>-a K-\a\\bA + T- 



1 + — a 
2' ' 



> 



-^{(H^,,,-2|^,|)(2 + aH) + 4r} 



4| l + 



2 

If (i?t t + a|/jtl) < 4a/?t,tT, then in this case: 



- (a^,,, H + /?,,,- a |Z7jf>0 (14) 



4| 1 + 



2 

Therefore, if |ci:|> 0 then for(/?t ^ + alZ^^ |) < 4a/?t ^r, we have A/> 0, the cost function increases. 
We now prove the Proposition 2: 

Proof: Let the ^th element 9^= 0 be activated as = 0 and the updated solution vector is denoted 
as X. The update of the cost function is then given by: 



A/ [x] -/ [x] = 1 9? - r ^ 



i+-kl (15) 

2l H 

Thus, if |xtp/?tt/2 + Refx^bt) ~ l^tl/(l + ^l^tl) < 0' then A/ < 0, there exists a nonzero 
value of the ^-th element that decreases the cost function that should be removed from the support. 

Table 2. HR-DCD Algorithm. 

Step Equation 

Initialization; x = 0, / = 0, r = y, b = A^r, R = A^A 

J Choose the first (^-th) element into the support according to: 

t = argmax;,|/?;,|2; / = {t],T = \bt;\ 

2 Repeat until a termination condition is met: 

3 Solve minx^ ||y - A/X/II2/2 + r^ne/ Unl/(1 + on the support /and 
update r using DCD iterations 

4 Update the regularization parameter: t <- yr, 0 < y < 1 

5 Remove the t-th element from / according to the rule: 

t = argmin;,e; \xj,\'^R],j,/2 + Reix^^bk] - t \xj,\/(il + a\X},\/2) s.t. 
\xt\^Rt,t/2 + Re{x;bt} - T |xt|/(l + a\xt\/2) < 0 
If the ^-th element is removed, update r: r r + x^A^^^ I <^ I/t 

6 Add the ^-th element into / according to the rule: 

t = argmaxk^ic (a\bk\ - Rkk)/aRl,, ■ [{Rk^k + albklY - 4aRk^kT^ s.t. (i?t,t + «^l^tl)^ > 4ai?t,tT 
If the t-th element is added, update: / <- / U t 



Combining the DCD algorithm in Table 1 and the above two propositions, we arrive at the HR-DCD 
algorithm presented in Table 2. 
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The HR-DCD algorithm starts with the zero support 1 = 0, and r = For each r, it updates the 
support and minimizes the cost function according to the corresponding propositions. Besides, the 
algorithm will terminate if the iteration times reach the preset parameter L or t < t^^^ (T^min is a 
predefined parameter). 

Remark 1: When the penalty is logarithmic function, we give two similar propositions as follows: 
Proposition 3: Add the ^-th (t G r) element into the support according to the rule: 

^ = argmax|Z7^| ^.f. |Z7j>r ^j^g^ 

Proposition 4: Remove the ^-th (t G /) element from the support according to the rule: 

^ = argimn^|x,|'/?,,+9i{x,\}-r-log(l + a|x,|) s.t. +9i{x>J-r-log(l + a|x,|) < 0 (17) 

Remark 2: When the penalty is 1^ norm, we can also give two similar propositions as follows: 
Proposition 5: Add the ^-th (t G /^) element into the support according to the rule: 

t = arg max \bj^ | s.t. Ib^f ^ > 2tR]~/ ( j g) 

Proposition 6: Remove the ^-th (t G /) element from the support according to the rule: 

r = argimn^|x,|'/^,,+9^{xX}-r|xJ sX.^\x]^ R^^^+^{x;b^]-T\xX <0 (19) 

The proofs for Equations (18)-(21) will be shown in the Appendix. By combining the corresponding 
propositions with related DCD algorithms, it is easy to obtain Hlog-DCD algorithm and H/^-DCD 
algorithm, respectively, and we omit the tabulated expressions here for simplicity. 

3.1. Complexity Analysis 

The complexity is given by P = QMN + 4NL + 2C^N + Q + 2^^^^!^ . The term QMN is for 
computing the initial b. The term ANL is for selection of elements in the support. The term 2Cy^N is for 
updating b in the total number of successful DCD iterations, and each update requires only 2N real 
valued additions, as multiplications by are bit-shifts. The term Q takes into account Q (total 
number of DCD iterations) tests to decide if the DCD iteration is successful. The debiasing (ZN^^ijLg) 
is now done by using extra Af^^^ DCD iterations on the finally identified support with size Lg. 

3.2. Off-Grid Problem Solution 

As stated before, we explore the idea of grid refinement technique [27] to reduce the number of the 
arbitrarily located potential positions of x. Two-dimensional passive radar imaging is selected as an 
example to illustrate the procedures of solution, which is realized as follows: 

(1) Under a coarse resolution Ax^, Ay^, where Ax^, Ay^ represent the range step and azimuth step 
respectively, calculate the coarse localizations (x^,y^),n = Ir" ,S results by finding 5 largest 
amplitude peaks using the proposed algorithms, i.e., HR-DCD, Hlog-DCD and H/p-DCD. 
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(2) Build a denser grid around tlie estimated locations (x^, with a finer resolution (Ax^, Ay^) as 

Figure 2 mentioned in Section 2. 

1 1 

(3) Set g = g + 1 and Ax^ = - Ax^, Ay^ = - Ay^, repeat for Q times from Step (1) until the desired 
resolution is achieved. 

In order to reduce the multiplications, we use dichotomy for grid refinement, which can use bit-shifts 
instead of multiplications to be implement in hardware system. According to our experimental results 
in Section 4, no more than 10-step grid refinement process is needed. 

3.3. Extension to Multiple Measurement Vectors (MMV) Case 

Although we have derived HR-DCD, Hlog-DCD and H/^-DCD algorithms from a single measurement 
vector (SMV), the results can be generalized to MMV case. In some kinds of applications, such as 
wideband source localization [42], multiple frequency bins are explored. As we all know that single 
frequency results can be easily disturbed by noise and interference, and in this issue, we can use a 
frequency diversity to achieve better performance. A traditional method of modifying the above 
approach to multi-frequency signals is averaging the results of all the frequencies. But due to the 
presence of noise, environmental mismatch and interference, different frequencies may give different 
results, thus this combining method may not work well. 

By the similar approach used in [28], our proposed methods require that for all frequencies to have 
the same support, which utilize the joint sparsity pattern, and choose an element to add to the support 
according to the corresponding propositions to achieve greater frequency diversity and avoid the possible 
wrong results. And the final result can be given by averaging all the results of all the frequencies through 
this way. 

4. Experimental Results 

In this section, we will present several simulation results which illustrate the effectiveness of 
the proposed methods. All experiments are performed by using MATLAB R2013a on a PC equipped 
with an Inter Core 17 3770k CPU, 3.5GHz and 32 GB memory. The state-of-the-art sparse recovery 
methods such as MP, BP, FOCUSS, SBL, H/l-DCD, etc. are selected for comparison. As a benchmark 
we will utilize an oracle sparse recovery (OSR) to represent the best inversion performance, which 
knows the true support of x. In addition, we set parameter a = 10 in Equation (4) for HR-DCD, in 
Equation (5) for Hlog-DCD, and p = 0.8 in Equation (6) for H/^-DCD, and = 10,Tn,in = 0.02, 
Y = 0.9 by experience from extensive simulation results. Other predefined parameters such as //, L, 5, 
and the grid refinement times Q, etc., are decided according to the actual conditions, which will be given 
in details in the next text. 

4.1. DO A Estimation 

As described by Malioutov et al, for one snapshot case, the DOA problem of Equation (5) in [27] is 
equivalent to Equation (2) in this paper. With the assumption that the sources can be viewed as point 
sources and their number is small, the underling spatial spectrum is sparse and it can be solved via 
sparse recovery methods mentioned above. The first simulation considers the scenario that there are 
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five uncorrelated signals impinging from [-47.7°, -28.5°, -9.2°, 11.6°, 30.1°]. The grid is set to be 
within the range of -90° to 90° with 1° spacing. A 30-element uniform linear array (ULA) spaced in 
half- wavelength units is used, and the signal-to-noise ratio (SNR) is set to be 20 dB. Here we only 
consider one snapshot case, and set // = 1,L = 5 = 8, Q = 10. Figure 4 depicts the solutions solved 
by different methods, and the vertical dashed lines mark the true directions. Obviously, our methods 
find the positions and the amplitudes exactly and achieve much sparser solutions than MP and BP. 
Moreover, they have a better amplitude estimation than H/l-DCD [40]. 

Figure 4. Solutions of MP (a), BP (b), H/i-DCD (c), HR-DCD (d), Hlog-DCD (e), 
H/p-DCD (f) for off-grid DO As. True DO As are denoted by the vertical dashed lines. 
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We give a time cost for the aforementioned methods in Table 3. From the results, we can see that 
our methods have the same time cost as the H/i-DCD. Although they are a little longer than MP, but 
they are much less costly than BP. However, the significant advantages of our methods are that they 
are very easy to implement in a hardware platform because only bit-shifts are needed instead of 
multiplications similar as [37^0], besides, they can achieve the similar reconstruction performance as 
BP without off-grid problem. 

Table 3. Time consumption of different methods. 



MP 


BP H/i-DCD 


HR-DCD 


Hlog-DCD 


H/p-DCD 


Time(s) 0.013 


0.721 0.130 


0.122 


0.129 


0.126 



Next we compare the mean squared error (MSB) of position and amplitude estimation performances 
by applying the proposed methods and several popular methods (MP, BP, FOCUSS, SBL, H/l-DCD, 
OSR), which are defined as follows: 



Position MSE = 



^position(x) - position (x 



||jf?05'/?/(9/z(x)|| 



(20) 



Amplitude _ MSE = - 



x-x 



(21) 



where x is an estimate of x. 

Figure 5 a compares the MSEs of DO A estimation results of different methods under varying SNR, 
which are averaged by 100 Monte Carlo trials. From the curves of MSE of position versus SNR shown 
in Figure 5a, it can be seen that H/i-DCD and the proposed methods have similar results in position 
estimation, and outperform their counterparts. However, our methods have better amplitude estimation 
accuracy than H/i-DCD from Figure 5b, and are even very close to the OSR performance. 

Figure 5. MSE of position v^". SNR (a) and amplitude v^". SNR (b). Uncorrelated source 
signals come from 47.7°, -28.5°, -9.2°, 11.6°, and 30.1°. Herein one snapshot is used. 
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The MSEs of DOA estimation versus number of signals shown in Figure 6 are obtained over 
100 Monte Carlo trials with SNR equals to 20 dB. We can clearly see that our methods perform best 
compared to other methods with the same parameters, especially on the amplitude estimation results. 

Figure 6. MSB of position v^". number of signals (a) and amplitude v^". number of signals (b). 
SNR is fixed at 20 dB. 
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4.2. Passive Radar Imaging 

In this part, we aim to test the performance of the proposed methods in a passive radar imaging 
system. Based on the final echo Equation (17) stated in [6] which is similar to Equation (2) in this 
paper, sparse recovery techniques can be utilized to solve the problem of passive radar imaging. In the 
simulation, we choose seven digital video broadcasting (DVB) satellites as opportunity transmitters and 
10 receivers, and the initialization grids are 21 X 41 represented for the imaging scene 20 m X 20 m 
(here we use the axis function in matlab to make the scene limited to be 12 m X 10 m so as to achieve 
a better visual effect). There are nine off-grid scatterers with different reflection coefficients in 
the scene of interest as the red circles illustrated in Figure 7. The number of frequency samplings 
corresponding to each transmitter is 10, and the SNR equals to 20 dB. Other simulation parameters are 
the same as displayed in [6]. Here we set // = 2, L = 5 = 10 and grid refinement times Q = 10. 

Figure 7a shows the poor imaging result by BP due to the reason that in practice the scatterers are 
not located exactly on the pre-discretized grids, and the echoes are mismatched with the measurement 
matrix, therefore the performances of CS -based reconstructions are generally unsatisfactory. Figure 7b 
is the imaging result by H/i-DCD, and we can see that it can seek the exact locations of the target, 
however the amplitude estimation of the reflection coefficients is not precise. 

Figure 7c-e represents the imaging results by the proposed methods. As expected, the proposed 
methods do not have the off-grid problem from the perspective of both the position and amplitude 
estimation results, which show the potential of our methods to be applied in practical passive 
radar system. 
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Figure 7. Passive radar imaging results by applying BP (a), H/i-DCD (b), HR-DCD (c), 
Hlog-DCD (d), H/p-DCD (e). Left and Right figures show the amplitude and position 
estimates, respectively. Red circles denote the true scatterers of the target, and blue crosses 
represent the estimated results. 
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Figure 8 shows the MSEs of position and amphtude estimation results versus SNR by applying 
different methods, which are averaged over 100 Monte Carlo trials. In the simulation, we have assumed 
that there are nine off-grid scatterers with unit reflection coefficients distributed randomly in imaging 
region. It is seen that our methods outperform the conventional methods (MP, BP, FOCUSS, SBL) in 
further improvement of MSEs of position and amplitude estimates when off-grid target emerges. In 
the meantime, they have better amplitude estimation result than H/i-DCD with the same parameters 
shown in Figure 8b. Moreover, the performances of our methods approach the OSR very closely under 
varying SNR, both in position and amplitude estimates. 



Figure 8. MSB of imaging results versus SNR. MSB of position v^'. SNR (a) and amplitude 
v^. SNR (b). 
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4.3. ISAR Imaging 

We utilize the CS ISAR imaging model as discussed in [43], and assume that the translational 
motion compensation has been already achieved and thus only the rotational motion is considered. 
The relation between the received signal and the complex reflection coefficients can be written as 
Equation (10) in [43], which is equivalent to Equation (2) in this paper. 

We use the quasi real data of an airplane (B-727) provided by the U.S. Naval Research Laboratory 
to test the feasibility and performance of the proposed methods, which is available on the website 
http://airborne.nrl.navy.mil/-'Vchen/tftsa.html. The stepped frequency radar operates at 9 GHz with the 
equivalent pulse repetition interval (PRI) of 3.2 ms and has a bandwidth of 150 MHz. For each pulse 
train, 64 complex range samples were saved, and the file contains 200 successive pulse trains within 
the long coherent processing interval (CPI). An additive noise is added to the original B-727s data, and 
the SNR is set to 10 dB. Herein we choose 40 range cell numbers and 20 cross range cell numbers 
for short CPI test, and the reconstruction of target spatial domain is discretized with 64 range bins and 
32 cross range bins, besides, we set // = 2,L = 5 = 60 and Q = S. Since in radar imaging field, the 
position estimates are very important when off-grid problem exists. 
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Figure 9 are the position recovery results by MP, SBL, FOCUSS, BP, Matrix-Pencil, ESPRIT, 
H/i-DCD and the proposed methods (HR-DCD, Hlog-DCD and H/^-DCD) shown as the red circles, 
and we use the conventional FFT-based ISAR image as the background for fair comparison. As stated 
above, the imaging performances of the classical methods (MP, BP, FOCUSS, SBL) are not good 
because of the off-grid errors. Matrix Pencil and ESPRIT can direct solve the positions of the plane, 
but their performances suffer from low SNR and few snapshots, therefore the imaging results are not 
good enough. In contrast, our methods can deal with the off-grid target, and achieve satisfying imaging 
performances which have better visual effect for the airplane shape. Moreover, it is obvious that the 
imaging results of our methods are better than H/i-DCD with the same parameters, especially on the 
details of airplane wings and tails as illustrated in zoomed-in regions of Figure 9h-j, which demonstrate 
the advantage of homotopy DCD method with non-convex regularization compared with convex 
regularization under the same measurements. 

Figure 9. ISAR imaging results of B-727 by MP (a), SBL(b), FOCUSS (c), BP (d), 
Matrix-Pencil (e), ESPRIT (f), H/i-DCD (g), HR-DCD (h), Hlog-DCD (i), H/^-DCD G). 
The red circles denote the corresponding position recovery result, and background represents 
the FFT-based reconstruction. 



10 



53 20 
E 



30 



S)40| 



50 



25' — 

^15 20 



60 



12 13 14 



10 



CD 20 I 

E 



30 



g)40| 



50 



r15 20 



60 



12 13 14 



10 15 20 25 
corss range cell number 

(a) MP 



30 



10 15 20 25 
corss range cell number 



30 



(b) SBL 



10 



o 20 1 

SI 



30 



.40 



50 



60 



25 ' — 

^15 20 



12 13 14 



50 [ I 
20 21 22 23 



10 15 20 25 
corss range cell number 

(c) FOCUSS 



30 



10 



CD 20 



30 



S)40| 



50 



60 



rl5 20 



12 13 14 



10 15 20 25 
corss range cell number 



30 



(d) BP 



Sensors 2014, 14 



5946 



Figure 9. Cont. 
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In order to quantitatively evaluate the amplitude estimation performances of the obtained ISAR 
images via different methods, we use the correlation value [44] to evaluate the similarity between the 
recovered images and the reference image, and the image entropy [44] to measure the focusing quality 
of the recovered images. They are defined as follows: 



£'nf(xJ = -^p,.log/?,. 



(22) 
(23) 
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where and denote the reference image vector using full measurements and the recovered image 
vector by different methods, and pi is the histogram of the recovered gray level (0^-255) image. 

The correlation and entropy values of the recovered ISAR images under different methods are 
summarized in Table 4. It can be seen that the proposed methods have a higher correction with the 
reference image and a lower entropy than their counterparts, and thereby exhibit a better capability of 
target-information extraction. 



Table 4. Correlation and entropy values by different methods. 





MP 


SBL 


FOCUSS 


BP 


Matrix 
Pencil 


ESPRIT 


H/i-DCD 


HR-DCD 


Hlog-DCD 


ffip-DCD 


Cor 


0.817 


0.890 


0.887 


0.913 


0.816 


0.798 


0.925 


0.952 


0.957 


0.965 


Ent 


0.826 


0.811 


0.836 


0.859 


0.810 


0.802 


0.713 


0.689 


0.683 


0.613 



5. Conclusions 

Two major problems of applying CS to sparse continuous signal reconstruction are discussed in 
this paper, and we have presented computationally efficient and accurate methods to overcome the 
difficulties. For solving the off-grid problem, we utilized the grid refinement technique. In the meantime, 
in order to obtain high recovery performance and keep low calculation complexity, we propose a fast 
and accurate homotopy DCD reconstruction combined with three typical non-convex regularizations, 
which promotes sparsity more strongly than any convex penalty function can. Extensive experiments 
have been conducted to validate and compare the performances of the proposed methods with 
several popular solvers. Our future work will try to synthesize with parallel sparse optimization 
technique using multi-core CPUs/GPUs [45], which may provide a viable solution to real-time 
potential applications. 
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Appendix 

Remark 1: the Penalty is Logarithmic Function. 

Proof for Proposition 3: Let the ^-th element = 0 be activated as = a: = |Q:|e^^^^^^\ and the 
updated solution vector is denoted as x. The update of the cost function is then given by 

A/ =/[x]-/[x] = i|6^p/^^^-9i|6^*^J + 2--log(l + a|6^|) (24) 

The cost function is reduced if A/ < 0. For a fixed A/ achieves a minimum if arg{a:) = arg {fo^}, 
and in this case 



A/ = -^\af' j — \a\\bj\ + T — log(l + a\a\^ 



(25) 



As we can see, if \a\ = 0, then A/ = 0. So if < 0, when \a\ > 0, there exists a \a\ that 

^'^1 |a:|=0 

makes A/ < 0. 

D I I IZ,I ^ 

dial l + a\a\ v^^; 



if^ 

d\a\ 



< 0, we can get 

|a|=0 



T-\b,\<0 (27) 
Therefore, if |£»tl > t, there exists a\a\ that decreases the cost function. 
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Proof for Proposition 4: Let the ^th element ^ 0 be activated as = 0, and the updated 
solution vector is denoted as x. The update of the cost function is then given by 

A/ = /[x]-/[x] = i|xj'i?,,,+9i{x>J-rilog(l + aK|) (28) 

Thus, if - |xtp/?t t + Refxt^t) ~ ^^'^^ ~^ ^l^tl) < 0' then AJ < 0, there exists a nonzero value 
of the ^-th element that decreases the cost function that should be removed from the support. 

Remark 2: The Penalty is Ip Norm Function. 

Proof for Proposition 5: Let the ^th element x^ = 0 be activated as x^ = a = |a:|e^^^^^^^, and the 
updated solution vector is denoted as x. The update of the cost function is then given by 

A/ = / [x] - / [x] = i \af R^ ^ - SR [ab^ } + ^ |^^r (29) 

The cost function is reduced if A/ < 0. For a fixed \a\, A] achieves a minimum if arg{a) = arg {b^}, 
and in this case 



1 . r u^w 

AT I |2 r» I Il7 I \ \P 

A/ = -cif K-\a\\h\ + T\a\ 



2' 



\a\ 



+ T\a 



2R. 



(30) 



R.. 



From Equation (30) we can see that if |a| = — and T\a\P - — < 0, then AJ < 0. 
Thus, the value of |a| resuhs in the decrement of the cost function is given by 

(31) 

And we can get the proposition that if l^t p"^ > 2tR]~^ , the ^-th element would be added into the support. 

Proof for Proposition 6: Let the ^-th element x^ ^ 0 be activated as x^ = 0 , and the updated 
solution vector is denoted as x. The update of the cost function is then given by 

Aj = j[x]-j[x] = ^\x^f R^^ + m {xA }-t\x^ r (32) 

Thus, if ^ |xtp/?t t + Re{x^bf^] — rlx^l^ < 0, then AJ < 0, there exists a nonzero value of the t-th 
element that decreases the cost function that should be removed from the support. 
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